*APCED - Analysis for community cohort's analysis
*Hyesung "Hace" Oh
*Started 11/07/2021
*Last worked on: 06/24/2022
****************************************************************************
*GENERAL*
clear all
set maxvar 7500

*DATE*
global currdate: display %td_CCYY_NN_DD date(c(current_date), "DMY")
global date = subinstr(trim("$currdate"), " ", "_", .)

*LOG*
cd "P:\apced\h2o\datastore\Aim2\source"
capture log close 
log using ..\output\analysis_comm_cohort_$date.log, text replace


	*Three-category : In addition to APC Involve
	*Low-volume (<8 E/M visits all_total)
	*High-volume, low APC (>= 8 visits all_total + share of APC < 50%)
	*High-volume, high APC (>= 8 visits all_total + share of APC >= 50%)
	*APCs are driving volume --> More E/M visits mean 
	*When there is high-volume, it does not matter whether it is high physician and high APC 
	*More APCs --> higher involvement of care 

*DATA*

*MAIN*
capture program drop main
program define main

	do ..\source\programs.do
	use final_comm_af.dta, clear
	alt_model_clean
	unadj_covar_assess
	gen_pinv_weights
	pinv_outcome_plots
	gen_pvol_weights
	areg_trt_models
	areg_alt_models //Run and store estimates
	output_results_1
	var_scaling
	ipwra_models //Run and store estimates -- NEED TO WORK ON THESE
	output_results_2 // Output areg and ipwra models
	other_sensitivity 
	output_results_3

	log close

end


****************************************************************************
capture program drop alt_model_clean
program define alt_model_clean

	capture drop high_volume
		gen high_volume = (all_total >= 8)
	
	capture drop apc_volume
		gen apc_volume = .
			replace apc_volume = 0 if high_volume == 0 
			replace apc_volume = 1 if apc_rate < 0.5 & high_volume == 1
			replace apc_volume = 2 if apc_rate >= 0.5 & high_volume == 1
	
	
	label variable apc_volume "High volume and APC rate"
	capture label drop a_v
		label define a_v 0 "Low volume (< 8 E+M visits)" 1 ">= 8 E+M visits& <50% APC rate" 2 ">= 8 E+M visits& >= 50% APC rate"
		label values apc_volume a_v
	tab apc_volume
	
	keep if apc_involve != . 
	
	*Generate an "all ER" variable
	capture drop all_er_count
		gen all_er_count = er_count + er_count_last30
		
	capture drop ever_base_er
		gen ever_base_er = er_count > 0
		tab ever_base_er
		
	capture drop er_base_decile
		xtile er_base_decile = all_er_count, nq(10)
		tab er_base_decile
		
	capture drop er_above_med
		sum er_count, detail
		local p75 = `r(p75)'
		gen er_above_p75 = (er_count >= `p75')
		tab er_above_p75
		
	*Generate percentage variable for outcomes: 30-day hosp and hospice death:
	capture drop hosptl_30day100
		gen hosptl_30day100 = hosptl_30day * 100
	capture drop hospice_death100
		gen hospice_death100 = hospice_death * 100
		
	
end

capture program drop unadj_covar_assess
program define unadj_covar_assess

	tab apc_involve

	foreach v in apc_involve apc_volume {
	    
	*Assess covariates
		table1by_gen age adrd_years female white black hispanic ccw dual_9mo high_volume all_total er_count acute_days ip_days ltc_days hh_days  hospc_days metro micro metro micro ahrf_rucc ahrf_totpop ahrf_pophisp10 ahrf_hcc pci medinc pov ahrf_hsdip ahrf_college zip_hmo_share zip_dual_share zip_sdi snf_bed hha hosptl_30day hosptl_death hospice_death using ..\output\table1_comm_by_`v', by(`v') mean_dec_place(2) sd_dec_place(2) replace
		
	}
	
end


capture program drop areg_trt_models
program define areg_trt_models

	tab high_volume apc_involve
	
	tab high_volume, sum(apc_rate)
	
	tab apc_involve, sum(high_volume)

	areg apc_rate high_volume i.year i.state_group age adrd_years female white black hispanic ccw dual_9mo high_volume all_total er_count acute_days ip_days ltc_days hh_days home_days hospc_days res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha, absorb(hkzip) cluster(hkzip)
	
end


capture program drop areg_alt_models
program define areg_alt_models

eststo clear

display "APC involve as treatment"
display "Then, APC/Volume as treatment"
display "First: full sample models"

	tab apc_involve

	foreach var of varlist hosptl_30day100 hospice_death100 {
	    
		eststo, prefix(`var'_areg_inv_): areg `var' er_count i.apc_involve i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha, absorb(hkzip) cluster(hkzip)
		
		eststo, prefix(`var'_areg_vol_): areg `var' er_count i.apc_volume i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha, absorb(hkzip) cluster(hkzip)
	
}	
	
end 


capture program drop output_results_1
program define output_results_1

	foreach var of varlist hosptl_30day100 hospice_death100 {
		
		esttab `var'_areg_* using "../output/`var'_all_areg.rtf", replace drop(*state_group* *year* er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days *res_rurality* ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha) b(%9.1fc) ci(%9.1fc) par scalars(F rmse) sfmt(%9.6fc) mtitles("APC Involve FS" "APC/Volume-Int FS")  label nogaps compress 
		
	}

end


capture program drop var_scaling
program define var_scaling


foreach v in er_count er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha {
    
	quietly sum `v'
	if `r(max)'>0 & `r(max)' !=. {
		replace `v' = `v'/`r(max)'
	}
	
	
}


end

capture program drop ipwra_models
program define ipwra_models

	eststo clear
	capture drop res_cnty_num
		destring res_cnty, gen(res_cnty_num)
	
	capture drop violator_*
	foreach var of varlist hosptl_cnt er_cnt {
		eststo, prefix(`var'_ipwra_): teffects ipwra (`var' i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha) (apc_involve i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha), osample(violator_inv_`var')
		
	*	drop if violator_inv_hosptl_30day100 == 1
	
		eststo, prefix(`var'_ipwra_): teffects ipwra (`var' i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha) (apc_volume i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha), osample(violator_vol_`var')
	
	}


end


capture program drop output_results_2
program define output_results_2

	foreach var of varlist hosptl_cnt er_cnt {
		
		esttab `var'_ipwra_* using "../output/`var'_IPWRAmodels.rtf", replace drop(*state_group* *year* er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days *res_rurality* ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha) b(%9.1fc) ci(%9.1fc) par scalars(F rmse) sfmt(%9.6fc) mtitles("`var': IPW RA w/ APC Involvement Trt" "`var': IPW RA w/ APC/Volume Trt") label nogaps compress
		
	}

end

capture program drop other_sensitivity
program define other_sensitivity 

	eststo clear

display "APC involve as treatment"
display "Then, APC/Volume as treatment"
display "Another alternative specification: High APC x High Volume interaction"

	capture drop apc_rate50
	gen apc_rate50 = apc_rate >= 0.5
	
	capture drop low_volume 
	gen low_volume = high_volume == 0
	
	capture drop highvol_apc50 
	gen highvol_apc50 = .
	replace highvol_apc50 = 0 if apc_rate50 == 0 & high_volume == 0
	replace highvol_apc50 = 1 if apc_rate50 == 1 & high_volume == 0
	replace highvol_apc50 = 2 if high_volume == 1
	
	capture drop lowvol_apc50_int
	gen lowvol_apc50_int =  low_volume * apc_rate50

	capture drop violator_*
	
		foreach var of varlist hosptl_30day100 hospice_death100 {

		eststo, prefix(`var'_areg_int): areg `var' er_count i.highvol_apc50 i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha, absorb(hkzip) cluster(hkzip)
		
		eststo, prefix(`var'_ipw_int): teffects ipwra (`var' i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha) (highvol_apc50 i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha), osample(violator_vol_`var')
	
	}

end


capture program drop output_results_3
program define output_results_3

	foreach var of varlist hosptl_30day100 hospice_death100 {
		
		esttab `var'_areg_int* using "../output/`var'_areg_int_models.rtf", replace drop(*state_group* *year* er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days *res_rurality* ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha) b(%9.1fc) ci(%9.1fc) par scalars(F rmse) sfmt(%9.6fc) mtitles("`var': AReg Int Trt") label nogaps compress
		
		esttab `var'_ipw_int* using "../output/`var'_ipw_int_models.rtf", replace drop(*state_group* *year* er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days *res_rurality* ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha) b(%9.1fc) ci(%9.1fc) par scalars(F rmse) sfmt(%9.6fc) mtitles("`var': IPWRA Int Trt") label nogaps compress
		
	}

end


capture program drop gen_pinv_weights
program define gen_pinv_weights


	capture drop pinv1
	capture drop pinv2
	capture drop pinv3
	mlogit apc_involve i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha, base(0) robust difficult	
	predict pinv*
	sum pinv1 pinv2 pinv3
	
	
graph box pinv1, over(apc_involve, label(labsize(small))) noout	///
title("{bf:Community cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene chars, state FEs", size(small)) title("{bf:Figure A1.} Predicted probabilities of being in the minimal APC involvement group")	///
ytitle("Predicted probability of being in" "{bf:minimal} APC involve treat group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(comminv_p1, replace)
	graph export ..\output\FigA1.tif, replace width(4000)

	
graph box pinv2, over(apc_involve, label(labsize(small))) noout	///
title("{bf:Community cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene chars, state FEs", size(small)) title("{bf:Figure A2.} Predicted probabilities of being in the moderate APC involvement group")	///
ytitle("Predicted probability of being in" "{bf:moderate} APC involve treat group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(comminv_p2, replace)
	graph export ..\output\FigA2.tif, replace width(4000)


graph box pinv3, over(apc_involve, label(labsize(small))) noout	///
title("{bf:Community cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene chars, state FEs", size(small)) title("{bf:Figure A3.} Predicted probabilities of being in the extensive APC involvement group")	///
ytitle("Predicted probability of being in" "{bf:extensive} APC involve treat group", size(medsmall)) ysize(9) xsize(16) ///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(comminv_p3, replace)
	graph export ..\output\FigA3.tif, replace width(4000)

end


capture program drop pinv_outcome_plots
program define pinv_outcome_plots

	capture drop pinv_tails
	sum pinv1, detail
	gen pinv_tails = (pinv1 < r(p5) | pinv1 > r(p95))
	
	graph twoway lpoly hosptl_30day pinv1 if apc_involve==0 & pinv_tails == 0 || lpoly hosptl_30day pinv1 if apc_involve==1 & pinv_tails == 0 || lpoly hosptl_30day pinv1 if apc_involve==2 & pinv_tails == 0, scheme(sj) legend(pos(1) ring(0) col(1) lab(1 "Low APC") lab(2 "Limited APC") lab(3 "Extensive APC")) xtitle("Predicted probability: Low APC") ytitle("Probability: Hospitalized during final 30 days") title("Pr(30-day hospitaliztion) vs. Predicted Pr(Low APC)")
	graph export ..\output\Comm_Hosp30_PropLowAPC.png, replace
	
	graph twoway lpoly hospice_death pinv1 if apc_involve==0 & pinv_tails == 0 || lpoly hospice_death pinv1 if apc_involve==1 & pinv_tails == 0 || lpoly hospice_death pinv1 if apc_involve==2 & pinv_tails == 0, scheme(sj) legend(pos(11) ring(0) col(1) lab(1 "Low APC") lab(2 "Limited APC") lab(3 "Extensive APC")) xtitle("Predicted probability: Low APC") ytitle("Probability: Death while enrolled in hospice") title("Pr(Death enrolled in hospice) vs. Predicted Pr(Low APC)")
	graph export ..\output\Comm_HospiceDeath_PropLowAPC.png, replace

end


capture program drop gen_pvol_weights
program define gen_pvol_weights


	capture drop pvol1
	capture drop pvol2
	capture drop pvol3
	mlogit apc_volume i.year i.state_group er_count acute_days age adrd_years female white black hispanic ccw dual_9mo ip_days ltc_days hh_days home_days hospc_days i.res_rurality ami anem asth atrfib chf chrkid copd depr diab glauc hipfx hyprlip bph hyperten hypothy ischhd osteo arthr stroke schi schiot bipl pvd ulcers psds drug cancer zip_sdi zip_hmo_share snf_bed hha, base(0) robust difficult	
	predict pvol*
	sum pvol1 pvol2 pvol3
	
	
graph box pvol1, over(apc_volume, label(labsize(vsmall))) noout	///
title("{bf:Community cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene chars, state FEs", size(small)) title("{bf:Figure A4.} Predicted probabilities of being in the low volume group")	///
ytitle("Predicted probability of being in" "{bf:low} APC vol group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(commvol_p1, replace)
	graph export ..\output\FigA4.tif, replace width(4000)

	
graph box pvol2, over(apc_volume, label(labsize(small))) noout	///
title("{bf:Community cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene chars, state FEs", size(small)) title("{bf:Figure A5.} Predicted probabilities of being in the high volume + <50% APC involvement group")	///
ytitle("Predicted probability of being in" "{bf:high} E+M <50% APC group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(commvol_p2, replace)
	graph export ..\output\FigA5.tif, replace width(4000)

	
graph box pvol3, over(apc_volume, label(labsize(small))) noout	///
title("{bf:Community cohort}", color(black) size(medsmall))	///
subtitle("Multinomial regression: bene chars, state FEs", size(small)) title("{bf:Figure A6.} Predicted probabilities of being in the high volume + >50% APC involvement group")	///
ytitle("Predicted probability of being in" "{bf:high} E+M 50%+ APC group", size(medsmall)) ysize(9) xsize(16)	///
ylabel(#10, angle(0) labsize(small))	///
ysc(titlegap(2))	///
box(1, color(dkgreen))	///
graphregion(color(white))	///
name(commvol_p3, replace)
	graph export ..\output\FigA6.tif, replace width(4000)


end



main
	